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SAMENVATTING  (ONGERUBRICEERD) 

De  3-dimensionale  eindige  differentie,  djddomein  (finite  difference  time  domain,  FD-TD)  methode  is  een 
numerieke  methode  om  elektromagnetische  penetratie-  en  verstrooingqnoblemen  op  te  lossen.  De  methode 
maakt  gebruik  van  een  eindige  differentie  lepresentatie  van  de  tijdafhankelijke  vergelijkingen  van  Maxwell. 
Het  te  nrodelleren  object  wordt  opgenomen  in  een  rooster  en  de  tijd  wordt  opgedeeld  in  discrete  intervallen. 
Door  toepassing  van  de  eindige  differentie  vCTgelijkingen  op  deze  discrete  ruimte  te  herhalen  voor  elke 
volgende  tijdstap  simuleren  we  de  voortplanting  en  verstrooing  van  golven.  In  dit  nq)port  worden  de  3- 
dimensionale  eindige  differentie  tijd  domein  methode  en  de  hiervoor  benodigde  algorithmen  bdumdeld. 
Resultaten  worden  gepresenteerd  voor  een  perfect  geleidende  plaat,  kubus  en  wig  en  een  dielektrisch 
gelaagde  bol.  De  berekende  resultaten  komen  goed  overeen  met  experimentele  en  exacte  thewedsche 
resultaten. 


TNO  report 


P!»ge 

4 


ABSTRACTT 

2 

SAMENVATTING 

3 

CONTENTS 

4 

1 

INTRODUCTION 

5 

2 

THEORY  AND  BASIC  ALGORITHMS 

6 

2.1 

Ideas  behind  the  FD-TD  method 

6 

2.2 

The  algorithms 

7 

2.2.1 

The  Yee  algorithm 

7 

2.2.2 

Lattice  regions  and  plane  wave  source  condition 

10 

2.2.3 

The  lattice  truncation  condition 

12 

2.2.4 

Sinusoidal  steady  state 

14 

2.2.5 

The  near-to-far  field  transformation 

15 

3 

FD-TD  COMPUTED  SCATTERING  AND  PENETRATION  PROBLEMS 

17 

3.1 

Radar  cross  section  of  a  square  metal  plate 

18 

3.2 

Radar  cross  section  of  a  rectangular  metal  cavity. 

20 

3.2.1 

Vertical  polarisation  of  the  incidmt  wave 

20 

3.2.2 

Horizontal  polarisation  of  the  incident  wave 

23 

3.3 

Radar  cross  section  of  a  metal  wedge 

25 

3.3.1 

Vertical  polarisation  of  incident  wave 

25 

3.3.2 

Horizontal  polarisation  of  die  incident  wave 

28 

3.4 

Field  penetration  in  a  lossy  inhomogeneous  sphen 

30 

4 

CONCLUSIONS 

33 

LITERATURE 


34 


TNO  report 


Fige 

5 


1  INTRODUCTION 

Electromagnetic  pmetratioa  and  scattering  problems  ate  usually  very  coiiq)licated.  Only  a  few  of 
them  can  be  solved  analytically.  For  the  remainder,  one  has  to  use  numerical  mediods.  Various 
methods  are  being  studied  at  the  moment  One  of  die  various  methods  under  consideration  is  the 
Finite-Difference,  Time-Domain  (FD-TD)  method. 

The  FD-TD  method  uses  a  finite  difference  representation  of  the  time  depoident  Maxwell 
equations.  The  object  of  interest  is  embedded  in  a  lattice,  and  time  is  divided  in  discrete  intervals. 
By  qiplying  the  finite  difference  equations  for  every  time  step,  we  simulate  the  prt^agation  and 
scattering  of  waves.  By  enforcing  a  plane  wave  at  t  -  0,  subsequent  time-stq>ping  wiU  result  in  a 
steady  state  which  can  be  monitored  and  used  to  calculate  the  object's  radar  cross  section. 

The  main  advantage  of  the  FD-TD  method  is  its  sinqilicity.  The  method  can  easily  be 
implemented  for  arbitrary  complex  objects,  because  it  is  possible  to  designate  every  cell  in  die 
lattice  its  own  material  constants. 

Chapter  two  deals  with  the  theory  and  the  used  algoridinis  for  the  3-dimensional  FD-TD  method. 
It  explains  the  complete  method  including  boundary  and  source  conditicHis.  In  chapter  diree  some 
results  are  presented  and  conqiared  with  other  mediods  and  measurements.  The  used  objects  are 
perfecdy  conducting  plates,  cubes  and  wedges  and  a  dielectric  layered  qibeie. 
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2  THEORY  AND  BASIC  ALGORITHMS 

2.1  Ideas  behind  the  FD-TD  method 

The  aim  of  the  FD-TD  mediod  is  to  model  the  propagatitn  of  an  electromagnetic  wave  into  a 
volume  of  space  containing  a  dielectric  or  conducting  object.  By  time-stq>ping,  Le.  iq>eatedly 
implementing  a  finite  differaice  analog  of  Maxwell's  time-dqiendent  curl  equations  at  each  cell 
of  the  grid,  the  incident  wave  is  tracked  as  it  first  propagates  to  die  object  and  dim  interacts  with 
it.  The  time-stepping  goes  on  until  the  sinusoidal  steady-state  is  achieved  at  all  points  in  the  grid. 

Time-stqiping  is  done  by  a  finite-difference  procedure  due  to  Yee  [1].  This  procedure  requires  the 
positioning  of  the  components  of  E  and  R  about  a  unit  cell  of  the  lattice  as  riiown  in  figure  1,  and 
the  evaluation  of  E  and  H  at  alternate  half-time  steps.  In  this  manner  a  second  order  accuracy  in 
the  pace  and  time  increments  is  attained  by  using  centred  finite-difference  expressions  for  both 
the  time  and  pace  derivatives. 


Fig.  1: 


POiitiopiiig  of  the  field  components  in  a  pace  cell 


TNO  report 


f 


Pkge 

7 


Beside  the  Yee  pfxx:edore  there  are  a  number  of  other  procedures  diat  need  to  be  developed.  The 
first  (me  is  the  lattice  truncatum  conditicm.  At  the  edges  of  die  lattice  the  centred  finite-difference 
expressions  can't  be  used  so  that  alternative  eiqiressions  fn  die  edges  have  to  be  developed.  These 
expressions  are  called  lattice  tnincati<m  (xmditions  or  absorbing  boundary  conditions.  These 
conditions  must  assure  that  the  outgoing  waves  pass  die  lattice  truncation  planes  widi  a  minimum 
of  reflections.  To  minimiTe  c<m^ter  storage  the  conditions  need  to  be  positirmed  as  close  to  the 
object  as  possible. 

The  second  procedure  that  is  needed  is  the  one  that  assures  that  there  is  a  plane  wave  travelling 
through  the  grid.  This  pnxxdure  is  caUed  the  source  conditicm.  This  source  must  not  be  the  cause 
of  any  extra  reflections  as  will  be  the  case  when  the  plane  wave  is  sinqily  imposed  cm  a  certain 
plane  in  the  grid. 

A  diird  {H-ocedure  is  needed  to  calculate  the  radar  cross  section.  Once  the  fields  in  the  grid  have 
reached  the  steady  state,  the  near  fields  are  available.  These  fields  will  have  to  be  used  to 
calculate  the  radar  cross  section.  This  includes  a  near-to-far  field  transformation.  In  the  following 
paragr^hs  a  detailed  discussion  is  given  of  the  various  algorithms. 

2.2  The  algorithms 


2.2. 1  The  Yee  algcsithm 

Assuming  that  e,  a,  p.  and  p'  are  isotropic,  Maxwell's  etpiations  in  Carthesian  (xxndinates  can  be 
denoted  as  follows 


1 

'dEy 

p 

dy 

as. 

1  dx  ■ 

dz 

-P’Hy] 


(la) 

(lb) 


(ld) 

(le) 

(lf) 


Where  e  is  the  electrical  permittivity  ia  faradsAneter,  o  is  the  electrical  conductivity  in 
mhos/meter;  p.  is  the  magnetic  permeability  in  henrys/meter;  and  p’  is  an  equivalent  mimetic 
resistivity  in  ohms/mete*.  The  magnetic  resistivity  term  is  provided  to  yield  symmetric  curl 
equations,  and  to  allow  for  the  possibility  of  a  magnetic  field  loss  mechanism.  In  accordance  with 
Yee  a  point  (x,y,z)  =  (iS  jSjkS)  in  ^ace  is  doioted  as  (iJJt)  in  which  6  s  5x  =  5y  s  Sz  is  the  space 
increment  of  a  cubic  lattice.  A  function  F(i5j5jr54)At)  of  space  and  tune  is  denoted  as  F°(iJ  Jt)  in 
which  At  is  the  time  increment  and  i,  j,  k,  and  n  are  integers.  Yee  used  centred  finite-difference 
expressions  for  space  and  time  derivatives  that  are  second  order  accurate  in  5  and  At  respectively. 


9F"(i.j,k)  ^  F"(i-K>^.j.k)-F"(i->^.j,k)  ^  2 

dx  8 

3F"(i.j.k)  ^  F"'*'>^(i.j,k)-F"->^(i.j,k)  ^  2. 

dt  At 


(2a) 

(2b) 


Using  these  expressions  to  get  a  finite-diffimnce  analog  of  the  MaxweU  equations  results  in 

Er‘(i  +  M.j.k)  =C3(m)E;(i+)^.j,k)  (3a) 

+C4(m)  (H^''^(i  +  ^,j  +  >^,k)-H“-'^(i  +  ^.j->^,k) 

-H;^^(i  +  j,k  +  +  K.  j.k  -  K)) 

E;i"‘‘(i.j  +  ^.k)  =C3(m)E5(i.j  +  >^,k)  (3b) 

+C4(m)  (H;-'^(i. j +y^,k-^y^)-  H;+^(i, j +y,k-y2) 

-H^*>^(i+ K.  3+ >^.k)+H^+>^  (i  -  3  +  y,k)) 


TNO  report 


Er'(i.j.k  +  >^)  =C3(m)E”(i,j.k  +  >^) 

+C4(m)  +  M. j,k  +  >^)  -  j.k  +  >^) 

j + K.k + Viy+nl*^  (i.  j  -  K.k + >^)) 
j  +  >^.k  +  >^)  =  Ci(m)H“‘-^(i.j  +  >^,k  +  ji) 

+C2(m)  (E“(i. j  +  3^,k  + 1) - Eja, j  +  )^,k) 

-Ej  (i.  j  +  kk  +  )^)+Ej  (i.  j,k  +  >^)) 
j.k  +  >^)  =  Ci(in)H;'>^(i  +  j.k  +  >^) 

+C2(m)  (E;(i  +  l,j,k  +  >^)-E“(i,j,k  +  )^) 

-Eja  +  j.k  +  1)+E“(i  +  K,  j.k)) 

Hf  ^(i  +  >^,  j  +  >^.k)  =  Ci(ni)H"-^(i  +  M.  j  +  K.k) 

+C2(m)  (E;(i  +  >^.j  +  l.k)-E“(i  +  >^.j,k) 

-Ej(i  + 1.  j  +  >^,k)+E^(i,  j  +  >^,k)) 


Ci(m)  = 


C2{in)  = 
C3(in)  = 


C4(m)  = 


2^.(m)- Aip'(m) 
2n(ni)  +  Atp'(ni) 
2At 

5(2^(ni)  +  Atp'(m)) 
2e(m)  -  Ata(in) 
2E(m)  +  Ato(in) 

2At 

5(2e(m)  +  AtO(m)) 


in  which  m  =  niedia(ij  J^)  is  die  type  of  medium  at  a  field  conqwnent  location.  To  ensure  stability 
of  this  time  stepping  algorithm.  At  is  chosen  to  satisfy  die  inequality  [2] 


where  c — is  the  maximum  wave  phase  velocity  widiin  diis  modeL  5  is  chosen  to  be  small 
compared  to  die  wavelength,  usually  5  S  VlO.  Also  5  has  to  be  small  conqiaied  to  the  dimensions 
of  the  object,  because  curved  boundaries  are  focmed  by  a  staircase  iqiproximation. 


TNO  report 


Page 

10 


With  the  system  of  equations  (3a)-(3f),  the  new  value  of  a  field  vector  conqwnent  at  any  lattice 
point  depends  only  on  its  previous  value  and  on  the  previous  values  of  the  conqionents  of  the 
other  field  vector  at  adjacent  points.  Therefore,  at  any  given  time  step,  the  computation  of  a  field 
vector  may  proceed  either  one  point  at  a  time,  or,  with  a  parallel  processing  conqnita',  at  many 
points  at  a  time. 


2.2.2  Lattice  regions  and  plane  wave  source  condition 

As  shown  in  figure  2  the  FD-TD  lattice  is  divided  into  two  distinct  regions,  separated  by  a 
rectangular  surface  u^iich  serves  to  connect  fields  in  each  region. 


Plan*  WiT«  Boundary 


AboerMnp  Boundary  Condttlon 


(A) 


Rg.  2:  Division  <rf  the  lathee  into  total>field  and  scattered-fiekl  regions  (a)  global  geranehy,  (b)  local  geometry 
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Region  1  of  the  lattice  is  dmoted  the  total-field  region,  region  2  die  scattered-field  region.  The 
object  is  posititmed  in  regitHi  1  where  all  field  quantities  are  conqirised  of  the  sum  of  die  iucident 
wave  and  the  scattered  field.  In  regicn  2  all  field  quantities  are  comprised  only  of  the  scattered 
field. 

The  rectangular  planes  that  fonn  die  boundary  between  regions  1  and  2  contain  Band  Hfield 
conqionents.  These  field  ccmqionents  are  updated  according  to  equations  (3)  and  subsequoidy 
corrected  to  maintain  the  two  distinct  regions.  Typical  FD-TD  confutations  at  diese  boundary 
points  are  as  follows. 

E;+‘(i. Jo.k  +  >^)  =  Ej+l(i. j„,k  +  K)|eqn(3c)  +C4(m)Hr^ (i. jo  "  +  K)  (6a) 

Er'(i  +  >^. Jo.k)  =  Er'(i  +  >2.  jo.k)|eqn(3.)  +  C4(m)Hr^ (i  +  >^, j„  -  ^.k)  (6b) 

+  >^)  =  -  >^.k  +>^)|«p(3e)  -C2(m)E‘"  (i,jo.k  +  ^)  (6c) 

+  ^2.  jo  -  >^.k)  =  (i  +  K. jo  -  VzMcqpOt)  - C2(®)Ex“  (i  +  K- jo.k)  (6d) 

Here  Ej^^HiJ^Jr+Vi)  is  the  usual  FD-TD  value  of  die  total  E,  cmtfonent  evaluated  at  point 
(iJoJr+Vi)  and  time  step  n+1.  The  superscr^  T  denotes  the  known  incident  field  component 
value.  These  computations  assure  ctmsistency  of  die  subtraction  pperatitms  of  field  components 
across  the  Region  1  /  Region  2  boundary.  In  effect,  total-field  quantities  are  always  subtracted 
fi-om  similar  total-field  quantities;  the  same  goes  for  scattered-field  quantities.  This  enforcement 
of  consistency  serves  to  precisely  connect  the  two  regions.  Further,  die  inclusion  of  arbitrary 
values  of  E*  and  H'  in  the  consistency  relations  permits  the  fecification  of  any  desired  plane 
wave  of  arbitrary  angle  of  incidence  and  arbitrary  polarisation. 

There  are  a  number  of  advantages  to  this  method.  For  the  complete  object,  the  high  dynamic 
range  total-field  formalism  is  retained.  This  permits  accurate  confutations  of  low-level  fields  in 
cavities  or  shadow  regions  of  the  scatteier.  A  very  accurate  simulation  of  die  radiation  condition 
is  possible  due  to  the  scattered-field  fwmalism  in  die  lattice  truncation  region.  An  otha  advantage 
is  that  the  incident  wave  contribution  need  to  be  confuted  only  for  die  field  confcments  at  die 
rectangular  surface  connecting  Regions  1  and  2.  And  last,  die  scattered  near  fieUi  in  Region  2  can 
be  easily  integrated  to  derive  die  far-field  scattering  and  radar  cross  section,  as  discussed  later. 
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2.2.3  The  lattice  ttuncation  condition 

Most  of  the  electromagnetic  field  problems  are  "open"  problems,  which  means  that  the  fields  are 
not  restricted  to  a  certain  area.  To  solve  these  problems  with  the  FD-TD  method  die  field 
conqiutation  zone  must  be  limited,  otherwise  an  imlimifeH  amount  of  data  will  have  to  be  stored. 
The  conqiutation  zone  must  be  large  enough  to  enclose  the  object  of  interest,  and  a  suitable 
boundary  condition  on  the  outer  perimeter  of  die  conqiutation  zone  must  be  used  to  simulate  the 
extension  of  the  computation  zone  to  infinity.  These  conditions  are  called  lattice  truncation 
conditions  or  absorbing  boundary  conditions.  Care  must  be  taken  because  these  conditions  must 
not  cause  spurious  reflections  of  outgoing  scattered  waves  and  must  not  cause  numerical 
instability. 


The  wave  equation  in  free  space  is 


a^u  a^u  a^u  i  a^u  „ 

ax^  dy  di}  c^  at^ 


with  U  =  Ejj,  Ey,  Ej,  H,,  Hy,  H^.  In  other  words 


LU  =  0 


with 

L  =  d2+d2+d2--Ld? 

c 


and 


L  can  be  written  as  L*  ■  L-  with 

c 

and 


/  f  \2 


(7) 


(8) 


(9) 


(10) 


(11) 


(12) 
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It  can  be  proven  [3]  that 

L“U  =  0  (13) 

applied  on  all  field  conq>onents  at  x  »  0  assures  that  all  outgoing  waves  are  absorbed.  At  x  »  h  a 
similar  expression  can  be  found: 

L^U  =  0  (14) 

For  the  numerical  implementation  of  these  conditions  a  Taylor  approximation  of  the  square  root 
in  the  L  operator  has  to  be  made. 

Vl-S^  =1-)^S^  (15) 

This  results  in  six  conditions,  one  for  each  face  of  the  boundary  surface. 


a^u  1  a^u  c  a^u  c  a^u  „ 
axat  c  at^  2  ay^  2  az^ 

(x  =  0) 

(16a) 

a^u^ia^u  ca^u 
axat  c  at^  2  ay^  2  az^ 

(x  =  a) 

(16b) 

a^u  1  a^u  c  a^u  c  a^u  „ 
ayat  c  at^  2  ax^  2  az^ 

(y=0) 

(16c) 

a^u  1  a^u  c  a^u  c  a^u  „ 
ayat  c  at^  2  ax^  2  az^ 

(y  =  b) 

(16d) 

a^u  1  a^u  c  a^u  c  a^u 

dzdt  c  dt^  2  dx^  2  dy^ 

(z  =  0) 

(16e) 

a^u  la^u  ca^u  ca^u  „ 
azat  c  at^  2  ax^  2  ay^ 

(z  =  c) 

(16f) 

The  numerical  implementation  of  these  conditions  further  involve  the  transfonnation  of  the 
differential  operators  to  finite  difference  expressions.  This  is  done  according  to  Mur  [4].  Typical 
transformations  are  (for  x  a  0); 


d^u"(K.j.»^) 

dxdt 


U"'^^(l.  j,k)  -  j,k)  -  U”~kl.  j.k)  +  U”-V0.  j.k) 

2-6At 


I 


r'  I 


r 


r 


1 


I 

I 


(17a) 
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a2u"(K.j,k) 

at2 


9V(K.j,k) 

ay2 


3^U"()^,j,k) 

dz^ 


U°'*'Ho.  j.k)  -  2U"(0.  j.k)+ U°-^(0,  j,k) 
2At^ 

.  u”-^^i.j.k)-2u°aj,k)-m°-^i.j.k) 

U°(0.  j  -H.k)  -  2U°(0.  j.k)  ■>•  U“(0.  j  -  l.k) 
25^ 

.  U“(l. j  +  l,k) -  2U“(l,i,k)  +  U”(l, j-  l.k) 
28' 

u°(0.  j.k + 1)  -  2U°(0.  j.k)-m“(0.  j,k  - 1) 
25^ 

.  U” (1. j,k  + 1)  -  2U“(1, j,k)  +  U“(l. j.k  - 1) 
2S2 


(17b) 


(17c) 


(17d) 


Finally,  the  absorbing  boundary  condition  at  x  =  0  becomes 

u"^‘(o.  j.k) = -u"-i(i,  j,k) + •  (u“+ki.  j,k) + u"-^o,  j.k)) 

cAt  +  S  '  > 

•(U“(0,j  +  l,k) 


(18) 


28 


—  •(U"(0.j.k)  +  U“(l,j.k))  +  : 


(cAt)" 


cAt  +  6  \  '  >  25(cAt  +  5) 

+U''(0,  j  -  l,k)  +  U"(l,  j  +  l,k)  +  U"(l.  j  -  l,k)  +  U“(0,  j.k  + 1)  +  U“(0,  j,k  - 1) 
+U"(1,  j,k  + 1)  +  U"(l,  j,k  - 1)  -  4  •  U"(0.  j,k)-  4  •  U"(l,  j,k)) 


For  the  other  absorbing  boundary  conditions  similar  expressirais  can  be  derived. 

2.2.4  Sinusoidal  steady  state 

Time  stepping  is  continued  until  the  desired  sinusoidal  steady  state  bdiaviour  is  achieved.  The 
time  needed  to  reach  this  state  is  mainly  depending  cm  the  object's  electrical  size.  The  total 
number  of  time  steps  needed  for  a  wave  travelling  at  die  qpeed  of  light  to  make  two  coaq>lete 
front-to-back-to-firont  traverses  is  normally  sufficimt  to  teach  die  steady  state.  Typical  times  can 
be  found  in  [S].  After  reaching  the  sinusoidal  steady  state  the  magnitude  ""H  phase  of  die 
components  of  the  near-fields  have  to  be  determined.  Due  to  die  fact  that  the  FD-TD  method 
causes  stmie  fields  to  have  a  noiqihysical  dc  offset,  an  algoridim  is  used  diat  determines  die 
maximum  and  minimum  of  a  field  coiiqioiieot.  The  rm^nitude  is  calculated  as  (max-min)/2  and 
the  phase  is  the  phase  of  the  maximum  minus  an  arbitrary  fixed  phase.  To  achieve  a  greater 
accuracy  multiple  cycles  can  be  treated  this  way  and  an  average  result  can  be  calculated. 


TNO  report 


Plage 

15 


2.2.5  The  near-to-far  field  transfonnation 

The  far  field  data  can  be  obtained,  in  princq>le,  by  solving  an  integral  equation  for  die  induced 
currents  on  the  surface  of  the  object.  For  complex  objects  this  can  be  very  difficult.  A  good 
alternative  is  to  set  up  an  equivalent  problem.  The  object  is  suironnded  widi  an  arbitrary,  closed, 
virtual  boundary,  the  surface  S, ,  on  which  the  near-field  data  is  obtained  via  the  FD-TD  method. 
This  surface  is  situated  in  the  scattered-field  region.  The  near-field  data  is  used  to  calculate 
equivalent  electric  and  magnetic  surface  currents  on  the  virtual  surface  S,. 

jseq(r)  =  “xH,(r)  (19a) 

Mscq(r)  =  -nxEs(r)  (19b) 

A 

where  n  is  the  outward,  unit  normal  vector  at  the  surface  S,,  see  figure  3.  By  making  the  interior 
of  surface  S,  enqity  with  zero  fields  and  no  sources  the  equivalent  problem  is  set  up.  The 
equivalent  surface  currents  on  produce  the  same  scattered  field  external  to  S,  as  in  the  original 
problem. 


(E  fi  ) 

region  B  •  • 


^  ) 

•  •  legion  B 


Fig.  3:  The  equivalent  problem 

The  scattered  far  fields  are  given  by  the  transform  of  the  equivalent  currents  of  (19a)  and  (19b) 
over  the  free  space  Green's  function  I6],r7].  If  (eo4i^  are  Ihe  region  B  medium  characteristics 
*  iJtA,  and  the  free  space  inqtedance  q^slZOlt  for  0  and  0  polarisations,  the  following 
scattered  far  fields  e^qnessions  are  obtained. 
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Ee  =  -jkollo 


E^=-jko^O 


Ae  + 


^0. 


A<^ 


^0 


where 

Aq  =  AxCos8cos<p+Ay  cos0sin<p-Aj8m9 
Fg  =  Fx  cos0cos<p+Fy  cos08in9-F2  sin0 
A^  =  -Ax  8in<j)+Ay  COS9 
F^  =  -Fj  sm9 + Fy  CO89 


and  the  potentials  in  the  far  field  region  are  given  by 


a' 

^Seq 

F 

4icr  JJ 

Mseq_ 

gjkorcos^^jg^ 


with 

rcos^  =  (xcos<ti  +  y  sin^)sin0  +  zcos0 


The  radar  cross  section  RCS  is  calculated  as  the  ratio 


RCS  =  47a2 


^_Ei  +  E|_^ 


r~4eo 


(20a) 

(20b) 


(21a) 

(21b) 

(21c) 

(21d) 


(22a) 


(22b) 


(23) 


where  E*e  and  E>^  are  the  corresponding  components  of  the  incident  plane  wave.  The  conqtlete 
bistatic  RCS  is  a  natural  result  of  this  procedure  for  a  given  incidrat  angle  and  polarisation  of  die 
illuminating  wave. 
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3  FD-TD  COMPUTED  SCATTERING  AND  PENETRATION  PROBLEMS 

The  used  3-diinensioiial  FD-TD  code  is  a  Fortran  code.  It  is  used  on  a  Convex  C230 
siqtercomputer  with  vectorisation  and  parallelisation  turned  on.  Typical  computation  times  for 
one  RCS  value  ate 

T  =  1.3  •  N  •  (1.7  •  10“®D^  +  2.7  •  lO'^D^  +2.9  •  10~^d)  (24) 

where  T  is  the  confutation  time  in  seconds,  N  is  die  total  number  of  time  steps  and  D  is  the  cube 
root  of  the  total  number  of  grid  cells.  The  memory  requiremmts  are  enf  irically  determined  at 

M  =  7.3  10~^D^  +  1.2  10"^D^+0.8  (25) 

where  M  is  the  number  of  megabytes  needed  by  the  FD-TD  code  for  single  precision  calculations. 
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3.1  Radar  cross  section  of  a  square  metal  plate 

In  diis  section  results  are  shown  for  the  FD-TD  con^nited  monostatic  radar  cross  section  of  a 
perfectly  conducting  square  plate.  The  results  are  coixq>ared  with  die  results  of  a  code  that  uses 
physical  optics  and  difhactitm  dieory. 

The  plate  spans  6A,  x  6X  and  as  8  is  chosen  VlO  the  plate  is  formed  hy  60  x  60  x  1  cells.  The 
absorbing  boundary  is  located  at  a  uniform  distance  of  10  cells  from  the  plate  surface.  The 
geometry  of  the  plate  is  shown  in  figure  4.  Verticai  polarisation  of  the  incident  wave  is  used. 
Azimuth  angles  ate  from  0*  to  10*  with  the  elevation  angle  fixed  at  90*. 


Y 


Fig.  4:  Oeometiy  <rfthe  plate 


The  results  ate  presented  in  figure  5,  whidi  shows  excellent  agreement  between  the  two  medKids 
at  the  maxima  within  1  dBsiiL  At  the  negative  going  peaks  a  greater  deviation  is  noticed,  due  to 
the  impel fecdons  embedded  in  the  two  methods.  However,  for  practical  applications  only  the 
maxima  and  the  positions  of  the  minima  are  important,  not  the  absolute  values  trf  die 


—  ph.  optics  &  edge  diff  •  FD-TD,  delta=lambd/10 
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3.2  Radar  cross  section  of  a  rectangular  meUl  cavity. 

3.2.1  Vertical  polarisation  of  tbe  incident  wave 

We  next  «ini»i<ter  die  target  shown  in  figure  6.  Tbe  shown  cavity  is  formed  by  a  perfectty 
conducting  cube  widi  (Hie  side  removed. 


Fig.  6:  Oeometty  of  the  metal  cavity 

The  cavity  is  ilhiminated  by  a  vertical  polatised  plane  wave  widi  a  fiequeocy  of  IS  GHz.  The 
dimensions  of  the  cavity  as  shown  in  figure  6  ate  a  ■  19.9  mm  and  d  ■  0.96  mm.  Using  6  * 
0.4975  mm  *  X/40  and  forming  die  cavity  by  a  >  40  cells  and  d  ■  2  ceDs  b  a  good  diactetisbion 
of  die  cavity.  The  FD-TD  coaqmtBd  monostadc  radar  cross  sectioa  versus  arimulh  an^  is 
compared  to  measwements  made  in  die  anedwic  dnunber  fimiltty  at  FEL-TNO.  Measniemeats 
are  made  accordhif  to  [8].  The  caUbcated  background  level  in  die  anechoic  room  was  bekm  - 
SS  dBsm.  The  elevation  angle  is  fixed  at  90*. 
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RCS  (dB  m2) 


azimuth  angle  (deg) 


•  FD-TD,  delta=0.995mm  —  measurement  data 


Fig.  7:  MeasiireiiientvenasIl>TD  code  for  a  inetal  canity.  votkalpotaiiMttianctfiiicMienlwsve 

Hgure  7  shows  a  teasonaUe  agreement  between  FD-TD  data  and  measurement  data,  widun  3  dB. 
However,  care  must  be  taken  because  die  measurement  data  is  not  as  e;q>ected.  At  an  angle  of 
180*,  the  rear  of  dw  caviar,  die  object  bdiaves  as  a  fiat  conducting  plate  which  has  a  well  defined 
monottatk  radar  cross  section  qiproximately  -23  dBsm  whidi  shows  more  agreement  widi  dw 
FD-TD  data  than  widi  dw  measurement  data,  ^ipareiidy ,  dw  only  deviation  is  an  ofiEwt  of  about  - 
2  dBsm  on  dw  measurement  data  relative  to  dw  FD-TD  results.  For  dw  remainder  dw  patterns 
seem  to  be  matching.  This  oiCMt  is  probably  canaed  by  a  calibration  enor  for  foe  reference  target 
used  during  dw  measmenwats.  For  dw  calibration  a  qfoete  was  used  afoich  had  an  RCS  value  of 
only  4  to  8  dB  above  dw  background  levd.  This  can  result  in  a  wrong  RCS  vahw  for  dw  reference 
object  whidi  oqdains  dw  obaerved  offeet 

FD-TD  calculations  uaiog  8  »  0.993  mm  —  A/20  and  fomifaig  foe  cavity  Iqr  a  »  20  cells  and  d  ■»  1 
cell,  show,  when  conyated  to  the  former  FD-TD  reauH,  diat  this  change  in  reaohitioo  givea  only  a 
di||n  dianfe  in  RCS  panem.  Figure  8  shows  font  dwosing  5  >  X/20  seems  to  be  sufficient  to 
model  foe  caviqr  widi  good  accuracy. 


15  30  45  60  75  90  105  120  135  150  165  180 


azimuth  angle  (deg) 


FD-TD,  deltas 0.4975mm  “-FD-TD,  delta=0.995mm 


The  infhieace  of  the  grid  resolntioa  on  die  FD-TD  lesuhs  fior  a  metal  cavity,  vatical 
polarised  incident  wave 
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3.2.2  H<»i2oatal  polarisatkm  of  die  incident  wave 

We  now  creiaiAa-  die  same  cavity  as  used  before  iUuminated  by  a  horizontal  polarised  plane 
wave.  Using  S  -  V40  figure  9  shows  that  die  FD-TD  con^nited  radar  cross  section  patton 
shows  good  with  the  measurement  data  excqrt  for  an  offset  cm  the  measuremoit  data 

of  qjproximately  -7  dBsm. 

RCS  (dB  m2) 


azimuth  angle  (deg) 

'  FD-TD,  delta=0.995mm  “  measurement  data 

Fig.9:  Measuremeot  versus  FD-TD  code  for  a  metal  cavity,  horizontal  polarisatioa  of  die  mcident 

wave 


Again  it  is  assumed  that  die  measurement  data  are  mconect  due  to  die  calibration  atot  mentumed 
in  the  pievkHis  section.  A  clue  diat  states  diis  ooochision  is  die  fact  diat  die  radar  cross  section  at 
an  azhmdh  angle  of  180*,  which  is  perpendicular  incMeoce,  should  be  indqiendent  of  die 
polarisation  of  the  incident  wave  because  at  fois  angle  die  cavity  is  symm^ric  for  rotatkms  90* 
along  the  x-axis.  The  FE>-TD  calculations  are  in  agreement  widi  dus  foct  but  die  measurements 
show  a  ducrepancy  of  afanort  3  dBsm  between  horizontal  and  vertical  polarisations.  Assuming 
die  FD-TD  calrolations  are  correct,  a  comparison  is  made  widi  H>-TD  calculatuns  niring  a  grid 
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3.3  Radar  cross  secdoa  of  a  metal  wedge 

3.3. 1  Vertical  polarisatitm  of  incident  wave 

Hie  next  analyzed  target  is  a  perfectly  conducting  wedge  widi  an  opening  angle  of  10’  as  shown 
in  figure  11.  At  a  fiequency  of  IS  GHz  side  a  is  dimensioned  to  be  SX.  Choosing  8  to  be  X/20  the 
length  a  is  formed  by  100  cells.  Again  the  boundary  is  only  10  cells  away  from  the  outer 
dimensions  of  the  wedge. 


Rg.  11;  Geomedy  of  the  wedge 

The  coitral  axis  of  the  wedge  is  positiooed  paraUel  to  die  x-axis.  The  computed  FIVTD  result  for 
a  veitkal  polarised  incident  wave  is  shown  in  figure  12  together  widi  the  experimental  results 
obtained  at  die  anechoic  dumber  facility  at  FEL-TNO.  The  elevation  angle  is  fixed  at  90*. 
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RCS  (dB  m2) 


azimuth  angle  (deg) 

"  FD-TD,  delta= lambda/20  —  measurement  data 


Fig.  12;  Measurement  versus  FD-TD  code  fOT  a  metal  wedge,  vertical  polarised  incident  wave 

A  good  agreement  is  obtained  for  the  main  lobe  at  85’  and  the  first  second  and  diird  side  lobe.  At 
the  angles  around  30*  and  140’  the  FD-TD  method  gives  wrtmg  radar  cross  section  values.  These 
are  probably  due  to  the  fact  that  die  two  planes  forming  the  wedge,  are  malring  an  angle  with  the 
grid  of  only  S’.  Such  a  small  angle  has  a  very  pOOT  discrete  representation.  To  show  this  effect  die 
wedge  is  positioned  as  in  figure  13.  The  same  experiment  is  repeated  and  die  results  are  shown  in 
figure  14.  An  imfxovemrat  can  be  noticed  at  azimuth  angles  around  30’  and  140’.  This 
improvement  will  become  even  more  obvious  in  the  next  paragrqrii  as  the  wedge  is  iiiiiininatwd 
by  a  horizontal  polarised  plane  wave. 


15  30  45  60  75  90  105  120  135  150  165  180 

_ azimuth  angle  (deg) _ 


•  FD-TD,  delta=iambcia/20  —  measurement  data 


Measurenjoit  vemu  FD>1D  code  for  a  metal  wedge  wiifa  an  inqxoved  poaitioaing  in  the 
grid,  veitical  potariaed  incident  wave 
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3.3.2  Horizontal  polarisation  of  &e  incident  wave 

As  in  the  previous  cluster  the  wedge  is  analyzed  with  the  3D  FD-TD  code  in  two  ways;  One  with 
the  central  axis  of  the  wedge  parallel  to  the  x-axis  and  die  second  (me  widi  one  plane  of  die  wedge 
parallel  to  the  x-axis  as  in  figure  13.  This  time  a  horizontal  polarisatioD  of  die  incident  wave  is 
used.  The  results  are  shown  respectively  in  figure  IS  and  figure  16.  Again  good  results  ate 
obtained  for  certain  angles.  When  analyzing  these  two  figures  it  becomes  very  obvious  that  die 
positioning  of  the  wedge  in  the  grid  has  a  major  effect  on  the  coopted  RCS  value  at  those  angles 
where  interference  of  the  reflection  of  two  planes  of  the  wedge  becomes  dominant. 

RCS  (dB  m2) 


azimuth  (deg) 


•  FD-TD,  delta= lambda/20  —  measurement  data 

Measurement  versus  FD-TD  code  for  a  metal  wedge,  horizontal  polarised  incident  wave 


Fig.  15: 


azimuth  angle  (deg) 


•  FD-TD,  delta*  lambda/20  ~  measurement  data 


Fig.  16: 


Measuremoit  versus  FD-TD  code  a  metal  wedige  with  an  improved  positkn  in  the  grid, 

horizontal  polarised  incidem  wave 
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3.4  Field  pmetraticm  in  a  lossy  inhnnogeneoas  sjdieie 

The  last  object  to  be  analyzed  is  a  lossy  inbcmiogaieous  dielectric  ^bere  niiicb  consists  of  a  core 
and  a  shell.  The  core,  marked  in  figure  17  with  1,  has  the  following  characteristics:  e^i  »  72, 
Cl  -  0.9  S/m  and  ai  s  0.08  m.  For  (he  shell,  maiked  with  2,  these  values  ate:  e  15, 

-  0.05  S/m  and  -  0<15  m. 


Hg.  17:  Geometry  (rf  the  lossy  inhomogeoeoas  qjben 

The  incident  plane  wave  is  travelling  in  die  negative  z-dnectioo  and  has  an  E  oonqioiieiit  in  die 
negative  x-directum  and  an  H  conqKmeiit  in  die  y-direction  as  in  figure  17.  The  centre  of  die 
^  This  electromagnetic  problem  has  an  analytical  sohitkm  diat  can 

be  obtained  with  the  Mie  series  ^].  The  of  die  field  components  of  die  electrical  field 

m  die  qiheie  ate  compoted  widi  die  FD-'TO  code  and  con^iaied  widi  die  anatytical  solution.  A 
frequency  of  100  MHz  is  used  and  die  ^here  is  diacretised  at  two  resolutions,  S  >  0.01  m  and 
6-0.02  m. 


TNO  report 


abs(Ez  O'/m 


X-Xc  (cm) 


FD-TD,  clelta=1cm  +  FD-TD,  delta=2cm  —exact  solution 


ig.  18:  Ez  versus  x  mside  a  lossy  mhcmiogeneous  spheic,  y  y^  z  s  z 

abs(Ex)  (V/m) 
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-  Y-Yc  (cm) - 

FD-TD.delta*1cm  +  FD-TD,  cieltas2cm  —exact  solution 


Ex  vems  y  kiridc  a  loaty  ininaioiaieoas  qtae,  X  a  x,,  z  ■  z. 


TNO  report 


Page 

32 


abs(Ex)  (V/m) 


X-Xc  (cm) 


•  FD-TD,  delta=1cm  +  FD-TD,  delta=2cm  —exact  solution 


Fig.  20:  Ex  versus  x  inside  a  lossy  inbomogeswous  sphere,  y  ■  z  > 


These  figures  show  a  good  agreement  between  the  FD-TD  cooqiuted  results  and  die  analytical 
solution.  It  becomes  clear  diat  at  sufficiendy  high  resolutions  FD-TD  is  abte  to  model  field 
penetration  in  lossy  dielectric  objects  widi  high  contrasts. 
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4  CONCLUSIONS 

The  S-dimeasioiial  FD-TD  medwd  is  a  useful  mediod  for  analyziiig  problems  of  scattering  and 
penetration  of  electromagnetic  waves.  modem  sqMrcomputers,  scattering  from  objects  with 

maiiiniim  dimensions  of  leas  than  10  wavekngdis  can  be  analyzed  wifh  dtis  m^uxl.  Dqwnding 
on  die  cnwylexity  of  die  alufie  and  die  materials  used  these  dimensions  may  become  smaller. 
Good  results  are  obtaioed  for  5  «  A/IO  to  5  »  X/40  depending  on  the  canplexity  of  die  object  of 
interest.  For  perfectly  conducting  objects  diat  fit  exactly  in  die  lattice,  5  =  X/20  seems  to  be 
sufficient  to  model  die  RCS  with  1  dB  accuracy.  For  other  objects  a  higher  discretisation  has  to  be 
used. 
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